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A Sum-of-Squares approach to the Stability and Control of 
Interconnected Systems using Vector Lyapunov Functions* 
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Abstract — Stability analysis tools are essential to understand¬ 
ing and controlling any engineering system. Recently, sum- 
of-squares (SOS) based methods have been used to compute 
Lyapunov based estimates for the region-of-attraction (ROA) 
of polynomial dynamical systems. But for a real-life large scale 
dynamical system this method becomes inapplicable because of 
growing computational burden. In such a case, it is important 
to develop a subsystem based stability analysis approach which 
is the focus of the work presented here. A parallel and scalable 
algorithm is used to infer stability of an interconnected system, 
with the help of the subsystem Lyapunov functions. Locally 
computable control laws are proposed to guarantee asymptotic 
stability under a given disturbance. 

I. INTRODUCTION 

Toward the end of the nineteenth century the Russian 
mathematician A. M. Lyapunov [1] has introduced a number 
of powerful tools for the stability analysis of nonlinear 
dynamical systems. His stability results have been later 
generalized by Barbashin, Rrasovskii, and LaSalle, while 
control system engineers have used Lyapunov’s methods to 
design stabilizing feedback controllers — see [2] and refer¬ 
ences therein. Nevertheless, constructing a system Lyapunov 
function is usually a difficult task, which becomes daunting 
when the size of the dynamical system increases. 

A more practical approach is to define the Lyapunov 
function of the interconnected system as some function of 
the subsystem Lyapunov functions. There are different func¬ 
tional forms for the Lyapunov function of the interconnected 
system, such as a scalar Lyapunov function expressed as 
a weighted sum of the subsystem Lyapunov functions, or 
applications of vector Lyapunov functions and comparison 
principles [3]—[6]. Formulations using vector Lyapunov func¬ 
tions [7], [8] are computationally attractive because of their 
parallel structure and scalability. In [4], it was shown that 
if the subsystem Lyapunov functions and the interactions 
satisfy certain conditions, then application of comparison 
equations [9]—[ 11] can provide a certificate of exponential 
stability of the interconnected systems. In this work we seek 
an algorithmic certification of asymptotic stability via the 
vector Lyapunov function approach, where each subsystem 
Lyapunov functions are expressed in some polynomial form. 


Primarily we will focus on an example of a randomly 
generated network of nine modifiecf] Van der Pol oscillators. 
Each Van der Pol oscillator can be represented as a two- 
state system with state dynamic equations as polynomials 
of degree three [12]. The network is then decomposed into 
many interacting subsystems. Each subsystem parameters 
are so chosen that individually each subsystem is stable, 
when the disturbances from neighbors are zero. Sum-of- 
squares based expanding interior algorithm [13], [14] is 
used to obtain estimate of region of attraction as sub- 
level sets of polynomial Lyapunov functions for each such 
subsystem. Finally a sum-of-squares based scalable and 
parallel algorithm is used to certify stability in the sense 
of Lyapunov of the interconnected system by using the 
subsystem Lyapunov functions computed in the previous 
step. A distributed control strategy is proposed that can 
guarantee asymptotic stability of the interconnected system 
under given disturbances. Following some brief background 
in Sec. [IT] we outline the problem statement in Sec. [Ill] 
An algorithmic approach to certifying asymptotic stability 
is presented in Sec. [Iv] while a distributed control strategy 
is discussed in Sec. El Sec. 1 VII shows an application of our 
stability analysis and control approach to a network of Van 
der Pol oscillators. We conclude the article in Sec. rvm 

II. BASIC CONCEPTS AND BACKGROUND 

Before formulating the problem, let us briefly review some 
of the key concepts behind our analysis. We will first discuss 
how the stability of a dynamical system can be analyzed by 
constructing suitable Lyapunov functions. Then we briefly 
refer to sum-of-square polynomials and a very useful result 
which helps us in formulating the sum-of-squares problems. 

A. Lyapunov Stability Methods 

Let us consider the dynamical system 

t >0, xeK", /(0, ix i) =0„ x i (1) 

We assume that x = 0, !X i (for simplicity, t will be often 
dropped when obvious) is an equilibrium point of the dy¬ 
namical systerr@, and / : R” —» R" is locally Lipschitz. An 
important notion of stability is as follows: 
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Definition 1 The equilibrium point at origin is called 
asymptotically stable if it is stable in the sense of Lyapunov, 

1 We choose the oscillator parameters in such a way that the Van der Pol 
oscillators have a stable equilibrium at origin. 

2 Note that this is not a restrictive assumption, since by shifting of state 
variables, the origin can always be made an equilibrium point. 


and if 


35 > 0 s.t. Ilx' 


<8 


; lim||x(f)|| 2 = 0. (2) 


The Lyapunov stability theorem [1], [15] presents a suf¬ 
ficient condition of stability through the construction of a 
certain positive definite function. 

Theorem 1 The equilbrium point x = 0„ x i of the dynamical 
system in 0 is stable in the sense of Lyapunov in 54 £ R", 
if there exists a continuously differentiable positive definite 
function V : 54 —> R (henceforth referred to as Lyapunov 
function) such that, 

V(0)=0 (3a) 

y(x) >0,Vxe^\{0„ xl } (3b) 

and, — V (x) > 0,Vx £ 54 (3c) 

Further, if V satisfies —V(x) > 0,Vx £ 54\{O kx i}, then the 
equilibrium point at origin is asymptotically stable in 54. 

When there exists such a function V (x), the region of 
attraction (ROA) of the stable equilibrium point at origin 
can be (conservatively) estimated as 

& A := {x £ 54 |y (x) < f nax } (4a) 

where, := argmaxjx £ R" |v(x) < 7 } C 54 (4b) 

It can be noted that, without any loss of generality, the 
Lyapunov function can be scaled by f nax , so that the ROA 
is given by, 

@ A :={x£R n \V(x)<l} (5a) 

where, V(x) = V(x)/f wx (5b) 

Henceforth, for simplicity, we would assume, without any 
serious loss of generality, that the ROA is estimated to be 
sub-level set of V (x) = 1. 

B. Sum-of-Squares and Putinar’s Positivestellensatz 

While the Theorem[I]gives a sufficient condition for stabil¬ 
ity, it is still not a trivial task to find a suitable function V (x) 
that satisfies the conditions of stability, even when the origin 
is a stable equilbrium point. Relatively recent studies have 
explored how sum-of-squares based optimization techniques 
can be utilized in finding Lyapunov functions by restricting 
the search space to sum-of-square polynomials [13], [14], 
[16], [17]. Let us denote as the set of all polynomials in 
x G R". Then, 

Definition 2 A multivariate polynomial p (x) EE S% n is a sum- 
of-squares (SOS) if there exist some polynomial functions 
hi(x),i = 1... r such that p(x) = 1 hf(x), and the set of 

all such SOS polynomials is denoted by 


E„ := {p (x) £ Sf. n | p is SOS } 


(6) 


toolbox SOSTOOLS [18], [19] along with a semidefinite 
programming solver such as SeDuMi [20]. 

SOS technique can be used to find a polynomial Lyapunov 
function V (x) : R" —» R, with V(0„ x i) = 0, which satisfies 
the following SOS conditions [13], [18], [19], [21]-[24], 


V(x)-0i(x) £ , Vx £ 54 

— y(x) — 0 2 (x) €E L„,Vx £ 54 


(7a) 

(7b) 


Given a polynomial p £ Sf n , checking if it is SOS is a semi¬ 
definite problem which can be solved with a MATLAB® 


for some domain 54 around the origin and positive definite 
functions 0 1 (x), 02 (x). Often it is convenient to choose 54 := 
{x £ R" |p(x) < ft }, for some positive definite function p(x) 
and p > 0. 

An important result from algebraic geometry called Puti¬ 
nar’s Positivstellensatz theorem [25], [26] helps in translating 
the SOS conditions into SOS feasibility problems. Before 
stating the theorem, let us define: 

Definition3 Given gj £ M n , for j = 1,2 ,...,m, the 
quadratic module generated by gj’s is (gi ,g 2 , ■ ■ ■ ,8m) '■= 
{ Oo + Tjj Gjgj | °o, Oj G Vj } 

Then the Putinar’s Positivestellensatz theorem states. 

Theorem 2 Let Jff = {x £ R" |gi (x) > 0,... ,g m (x) > 0 } be 
a compact set. Suppose there exists u(x) £ such that 

u(x)£Jt({gi,g 2 ,...,gm), (8a) 

and, {x £ R" |(/(x) > 0 } is compact. (8b) 

If p(x) is positive on Jd, then p(x) £ .M(g\.gi,... ,g m )- 

Note: often in this work, for the gf s used, the constraints 
® would be redundant, i.e. the existence of u(x) would be 
guaranteed [26fl 

III. PROBLEM OUTLINE 
Let us assume that the dynamical system in ( TJ) is in 
polynomial form, i.e. / is a vector of n polynomial^ Given 
the full dynamical system 0, we can decompose it into 
subsystems 

V/ = 1,2,.. .,m, 

Xi= fi(xi)+gi(x), Xj £ R"‘, (9a) 

/i'(0n,-xl) = 0„ ( xl, 5i(0nxl) = 0„ jX i (9b) 

X= (x\,xl,...,x T m ) T er (9c) 

m 

n = ^ rit , XiP\Xj = 0 (9d) 

i=l 

In the decomposed system description, the ffs denote the 
isolated subsystem dynamics, and gf s are the interactions 
from the neighbors [14], Let us also denote by 

-M := {«'} U {j | 3 xpXj, s.t. gij (xt,xj) / 0 } , Vi, (10) 

3 Simplicity in formulating an SOS problem motivated us to use Putinar’s 
version of Positivstellensatz over other, more general, versions of the 
Positivstellensatz theorem [26]. 

4 If the dynamics is not in polynomial form, it has to be recasted into a 
polynomial form, with possible additions of equality constraints [14], [21], 
[23], [27], a case not considered in this work. 



the set of neighbors (including the subsystem itself). We 
assume that the isolated subsystems are individually (locally) 
stable, and there exist Lyapunov functions for each of the iso¬ 
lated subsystems. The goal is to develop a framework for the 
stability analysis of the full interconnected system by using 
the local subsystem Lyapunov functions and considering the 
neighbor interactions. 

Given a decomposition <[9|, the next step is to find poly¬ 
nomial Lyapunov functions V}(x,): R™‘ —»R, such that [13], 

Vx ( - G S>i, Vi(xi) - 0n(x/) G (11a) 

and, - VV i (x,) r /}(x i -)-<fe(x,) g£„,. (lib) 
where, % := {x, G R" 1 ' |p,(x,) < ]3,} 

for some positive definite functions <j)n (x,), <j>a(xi)> p,(x,) 
and positive scalars ji,. Starting from an initial Lyapunov 
function candidate obtained using (ITU and a corresponding 
estimate of the region of attraction, an iterative process called 
expanding interior algorithm , [13], [14], is used to iteratively 
enlarge the estimate of the region of attraction by finding a 
better Lyapunov function at each step of the algorithm. At 
the completion of this iterative step, the stability of each 
isolated subsystem (assuming no interaction) is quantified 
by its Lyapunov function V,(x,), with an estimation of the 
boundary of the domain of attraction given by Si a a = 
{xi G R“<1 Vi( Xi ) < 1}. 

The Lyapunov level-sets can be used to express the 
strength of a disturbance. The equilibrium point of the system 
at origin corresponds to the level set V,(0„ ;X i) = 0, Vi. If there 
is a disturbance from this equilibrium point, the states of the 
system would move to some point x(0) away from the origin. 
This disturbed initial condition would result in positive level- 
sets V/(x,-(0)) = yP G (0,1] for some or all of the subsystems. 
A necessary and sufficient condition of asymptotic stability 
can then be translated into the condition 

Vi, V/(x,-(0)) = $ =► Vi, ; fim^(x;(f)) = 0 (12) 

In the rest of the article, we present SOS algorithms to 
test stability conditions and design local (subsystem-level) 
control laws to achieve asymptotic stability. 

IV. STABILITY UNDER INTERACTIONS 

It is assumed that the isolated subsystems in ([T]i are all 
(locally) asymptotically stable, and there exist subsystem 
Lyapunov functions V}(x,),Vt. The estimated region of attrac¬ 
tion of the interconnected system under no interaction, S? A , is 
given by the cross-product of the regions of attraction of the 
isolated subsystems, S#A.p which are defined as sub-unity- 
level sets of the corresponding (properly scaled) subsystem 
Lyapunov functions (as in (Q}-([5 ]i), i.e. 

f":=fA.iX%x-xf A ,»i (13a) 
where, S%A,i ={x; G R”' |V}(x,) < 1}, Mi (13b) 

In presence of non-zero interactions, the resulting ROA 
would be different. If there exists a Lyapunov function for 
the interconnected system, the ROA for the whole system 


could be expressed as some sub-level set of that Lyapunov 
function. While it is very hard to obtain a scalar Lyapunov 
function for the full interconnected system, one could use 
vector Lyapunov function approach to obtain certification of 
stability in a scalable way. 

In this present work, we choose not to impose any further 
restriction on the Lyapunov functions V)(x,-) than requiring 
that those are in polynomial forms, and concern ourselves 
with asymptotic stability. In general, it is difficult to test 
a necessary and sufficient condition of asymptotic stability, 
such as the one given in (IT4l) . Nevertheless, it is possible 
to derive sufficient conditions of asymptotic stability under 
certain scenarios. Let us now present a distributed iterative 
procedure which can be used to certify asymptotic stability 
in a domain defined by sub-level sets of the subsystem 
Lyapunov functions. 

A. Algorithmic Test of Asymptotic Stabiltiy 

Before proceeding to explaining our algorithm, let us first 
note the following result: 


Lemma 1 Suppose, for all i G {1.2,...,/;;}, there exists 
a strictly monotonically decreasing sequence of scalars 
{ef } ,k G {0,1,2,...}, such that 


Mi,k, Vj(x) := VV,(x;) r (/}(x,) +g,(x)) <0, Vx G @ k (14a) 


where, \= 



ef +1 < Vi{ Xi ) < ef, 
Vj(xj) < £y Mj ± i 


(14b) 


Then the system m is asymptotically stable in the domain 

{x G R" | nr=t Vi{xi) < £° }, if lim^ +00 ef = 0, vfl 


Proof: Please refer to Appendix [A] ■ 

Using the LemmaQ]we can devise a simple iterative SOS 
algorithm to certify whether or not a domain S> defined by 

Vi(xi') < v oi |, (15) 

for some scalars Voi £ (0,1],V/, is a region of asymptotic 
stability for the system in (0}. It is to be noted that, using 
the Putinar’s Positivstellensatz theorem (Theorem [2j, the 
condition in (fl4l) essentially translates into equivalent SOS 
feasibility conditions 



Mi,k, 3(71), ofj G En,-, s.t. 

- VV, r (fi + gi ) - 4 (Vi - ef +1 ) - £ 4 [e k j - Vj) 

where, «,• = V nj. 


06) 

G L,7 / 


The algorithmic steps to ascertain asymptotic stability are as 
outlined below: 

1) We initialize ef 1 = v a i,Mi G {1,2,... ,m}, and choose a 
sufficiently small £ G R + . 

2) At the start of the k -th iteration loop, we assume to 
know the scalars {e^e/,. ..,£*} ,V/, and our aim is 


5 If the limit condition does not hold, we can only guarantee stability in 
the sense of Lyapunov [1], [15] 




to compute the scalars ef +1 , Vi such that (IT6l > holds. 
Essentially we want to solve the optimization problem, 

Vi, min ef +1 (17a) 

rrk n -k c*+! 

°i0 ’ a ij’ b i 

s.t., the condition in (QU holds. (17b) 

This is solved by performing a bisection search for 
minimum e k+1 over the range [0, ef]. 

If ( fl7l) is infeasible at 0-th iteratiorfl for any i G 
{l,2,...,w}, we conclude that the system cannot be 
guaranteed to be asymptotically stable in S>, and abort 
the iteration. Otherwise we move on to step [3 

3) If (e* — ef +1 ) > £.Vi, we continue from step [2] for the 
(A+1 )-th iteration loop. Otherwise we stop the iteration 
deciding that the limits of the sequences {ef } , Vi, have 
been attained. Further, if the limits are all zero, we 
certify asymptotic stable in 3>. 

B. Remarks 

The algorithm presented in Sec. IIV-AI describes how one 
can determine asymptotic stability of an interconnected sys¬ 
tem in a domain 2> defined by the subsystem sub-level sets. 
This test can be performed locally, and in a parallel way, at 
each subsystem level. The Lyapunov functions Vi’s are to be 
found before the start of the analysis, and communicated to 
the neighboring subsystems. Then during each analysis, it is 
assumed that the neighboring subsystems can communicate 
with each other the computed sequences {ef} in real-time. 
With the help of the stored Lyapunov functions, and the 
updated {ef } of the neighbors, each subsystem will continue 
the iterative process outlined in Sec. IIV-AI Since only the 
neighbor information is required, this algorithm is reasonably 
scalable with respect to the size of the full interconnected 
system. Moreover, the algorithm motivates the design of 
a distributed control strategy that can ascertain asymptotic 
stability. 

V. DECENTRALIZED CONTROL 

In this section, we discuss design of a local and minimal 
control strategy such that the system in Q is asymptotically 
stable in a domain S> defined in ( IT5] >. We use the term 
minimal to suggest that the control be applied only in certain 
regions, and not everywhere, in the state space, while by the 
term local we suggest that the control be computable and 
implementable on a subsystem level. 

We envision the control to be computed by each subsystem 
at each iteration loop. At k-th iteration, Vk G {0,1,2,...}, 
the i-th subsystem, Vi G {1,2,... ,m} performs the following 
tasks: 

1) It identifies if it belongs to the following set 



r 

(fi + gi) >0, ) 

II 

i G {!,..., m} 

V- = £ k 

Vj<£ k '\/je^\{i} j 


(18) 

6 Because V/’s, ff s and g;’s are polynomial, if G) is feasible at 0-th 
iteration for all i, then it is also feasible for all subsequent iterations. 


which can be checked locally. If i £ f/ k , control is not 
necessary and it sets F k = 0„ | X i and proceeds to task [3] 
If, however, i G ff k , it proceeds to task [2] to compute 
a control law. 

2) If i G f/, a polynomial state-feedback control law F k : 
K"’ —>• R”‘, F k ( 0 njX i) = 0„ ;X i, is computed such that 



{v f =ef, Vj<e) V/e^U'}} 


<0. (19) 


This produces the equivalent SOS condition, 

-vy r (/i+g,+7* j-pf(ef-vi)- £ o£.(e?-^GZ* 

Pi e @m, °ij e l (1i Vj ^ i, n t = £ nj (20) 


3) Finally it performs the search over minimum ef +1 , as 
in 03, with the un-controlled subsystem dynamics 
(fi+gi) in the feasibility condition (IT6l) is replaced 
by the controlled dynamics (/) + gi + F k ). 

To summarize, each subsystem i computes control laws 
F k : R"‘ —> R"‘, with F k (0„. x i) = 0„. x i, during each k- 
th iteration, so that the subsystem dynamics under control 
becomes: 


Vi, VJt, Vx G , 


\ fi(xi)+gi(x)+F I k (x i ), 


where were defined in (1 1 4bb . 


i <t W k 

i G W k 


( 21 ) 


A. Remarks 

Often it is important to impose certain additional con¬ 
straints on the possible control laws, such as bounds on 
the control effort. Although control bounds can be easily 
incorporated in the SOS formulation, we decide to keep that 
for future studies. We note, however, that since we apply 
controls F k only on certain subsystems i G 6 // k , and in certain 
domains S> k C the control effort would be reasonably 
bounded. 


VI. RESULTS 

Let us describe the model of the interconnected system that 
we use here, and two examples to illustrate the applications 
of the stability analysis algorithm and control design. 

A. Model Description 

We will consider a network of nine Van der Pol oscillators 
[12], as shown in Fig. Q] The dynamics of each oscillator, 
in presence of neighbor interactions, is represented by 

Vy G {1,2, — ,9}, 

Xji = Xj 2 (22a) 

Xj 2 = PjXji (1 -xji) Xji +xj i £ £ jk x k2 (22b) 

kAi 

where jij’s are chosen randomly from (—2, 0) and the coeffi¬ 
cients, Cjk, of the interaction terms are chosen randomly from 
(—0.2, 0.2). It is to be noted, that the interactions need not 










together oscillators {2,3} and {5,6}, as shown below 



Fig. 1. A network of nine Van der Pol oscillators along with their isolated 
regions of attraction. 



Fig. 2. Comparison of estimated ROAs with true ROA for oscillator 1. 


be symmetric, i.e. in general, ^ Additionally, £,/. = 0 
if oscillator k is not a neighbor of oscillator j. Such choice 
of Hj’s ensure that the oscillators themselves are stable, 
with corresponding regions of attraction as shown in Fig. Q] 
The regions drawn in ‘red’ around each oscillator shows its 
true ROA, while the region in ‘blue’ shows an estimate of 
the ROA as sub-unity-level set of its polynomial Lyapunov 
function, as in (0. In Fig. [2] we compare with the true ROA 
the estimates obtained using a quartic Lyapunov function 
and a quadratic one. Also a sequence of estimates using 
the quadratic Lyapunov function are shown in ‘dotted black’ 
lines, which show how the ‘expanding interior’ algorithm 
iteratively expands the estimate of the ROA. Clearly, the final 
estimate improves as the degree of the polynomial Lyapunov 
function increases, but for this work we choose to stick to 
quadratic Lyapunov functions. 

We decompose the system into 7 subsystem^], by grouping 


7 This decomposition is arbitrary. [28] presents a method of decomposi¬ 
tion in weakly interacting subsystem, which however requires symmetric 
interactions. 


Si 

s 2 

S 3 

s 4 

s 5 

s 6 

s 7 


{osc 

1}; 


■A 

{S 2 ,S5,S 7 } 

{osc 

2, osc 

3}; 

^ 2 

{Si,53,56} 

{osc 

4}; 


^3 

{5 2 ,5 4 } 

{osc 

5, osc 

6}; 


{S 3 ,S 7 } 

{osc 

7}; 


^5 

{5i,56,5 7 } 

{osc 

8}; 


A 

{5 2 ,5 s } 

{osc 

9}; 


^7 

{Si,S 4 ,S5} 


( 23 ) 


Then we can write the subsystem dynamics, along with the 
neighbor interactions, in the form of ©. As an example, the 
states and the dynamics of S 2 are shown below, 

•*2 = / 2 (* 2 ) +g 2 (*); X 2 = (*21, *22, *31, X 32 f , 

*22 


h = 


82 = 


0 . 12 x 2 i* 3 2 -* 2 i — 0.4Lr 22 (l -x 71 ) 

X 3 2 

0.04x 3 lX 22 -X31 - 1.44x32 (l -Ajj) 
0 

—0.07xi 2 x 2 i 

0 

O.OIX12X31 +0.06x'4 2 X3i +0.1 x82X'31 


(24) 


Any randomly picked initial condition, x(0) € where 
is defined in d, can be mapped into corresponding 
subsystem Lyapunov function level sets, yf = V;(x,(0)),Vi. 
Then, by choosing v oi - = yf, we apply the iterative stability 
analysis algorithm to determine whether or not the domain 
in (fl5l > is a region of asymptotic stability, and if not, compute 
the necessary control by ( fl9l >. 


B. Example: Certifiably Stable without Control 

In Fig. |3(a)| the evolutions of all the states to an asymptot¬ 
ically stable initial condition is shown. The subsystem Lya¬ 
punov functions, shown in Fig. |3(b)[ monotonically decrease 
to zero starting from the initial level sets: 

7? = 0.242, y§ = 0.728, yf = 0.184, ^ = 0.658, 

7^ = 0.357, y£ = 0.279, y?= 0.283 (25) 

Then setting v oi - = yf, and choosing e = 0.001, we run the 
iterative stability algorithm which produces the results in 
Table J] At the end of the 2nd iteration, all the e}’s are zero, 
which certifies asymptotic stability of the full interconnected 
system. 


TABLE I 

Iteration results for a certifiably stable case 


k 

£ k 

£ k 


£ k 

b 4 


e 6 

p 

/ 

0 

0.242 

0.728 

0.184 

0.658 

0.357 

0.279 

0.283 

l 

0.024 

0.031 

0.147 

0.188 

0.004 

0.010 

0.127 

2 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 




























Fig. 3. Evolution of states and the subsystem Lyapunov functions for a certifiably stable initial condition. 




Fig. 4. Evolution of states and the subsystem Lyapunov functions for a non-certifiably-stable initial condition. 


C. Example: Certifiably Stable under Control 

Let us now present one example where the iterative 
algorithm fails to guarantee stability, and control is applied. 
Fig. [4] shows a stable initial condition, but the algorithm fails 
to certify stability for it. In Fig. |4(b)| Fi(/) is seen to increase 
initially before starting to monotonically decrease. When we 
apply the stability analysis algorithm to the initial level sets: 

y? = 0.953, = 0.990, $ = 0.149, $ = 0.479, 

^5 = 0.697, yg = 0.220, y? =0.103 (26) 

the iteration fails at the first iteration because the algorithm 
cannot find feasible £■ ’s for i = 1,2,5 (as shown in Table HU). 


TABLE II 

Iteration fails to certify stability 


k 

ef 

p 

p 

P 

t 4 

P 

p 

P 

0 

0.953 

0.990 

0.149 

0.479 

0.697 

0.220 

0.103 

1 

X 

X 

0.1010 

0.0056 

X 

0.0198 

0.0380 


Consequently decentralized control is activated for sub¬ 
systems Si,Si and S 5 . This results in certifiable asymptotic 
stability, with new ef's shown in Table [III] where the **’ 
denotes presence of controller^] which are documented in 

(ED. 

TABLE III 

Iteration certifies stability under control 


k 

P 

£? 

e| 

P 

t 4 

P 

P 

O 

P 

0 

0.953 

0.990 

0.149 

0.479 

0.697 

0.220 

0.103 

1 

0.004* 

0.004* 

0.1010 

0.0056 

0.003* 

0.0198 

0.0380 

2 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 


F 1 1 =[0,-0.569x 11 -2.271xi 2 ] r , 

7*2 = [0, —1 .237x22—0. 149x21, 0, —0.285x31—1 . 368 x 32 ] T , 
Fj= [0,-0.504x71 - 1.539x 72 ] r (27) 

8 In this example, we chose to seek linear controllers which proved to be 
sufficient. 




















































Fig. 5. Evolution of states and the subsystem Lyapunov functions, under decentralized control applied at Si ,£2 and S 5 , that certifies stability. 


It is to be noted that we decide to apply control only on the 
dynamics equations of the states xi 2 ,X 22,*32 and X 72 . Fig. 0 
shows that under the action of the controllers in (f27l) . all 
the subsystem Lyapunov functions decrease monotonically 
to zero. 


VII. CONCLUSIONS 

In this work, we present an algorithmic approach to cer¬ 
tify asymptotic stability of an interconnected system whose 
dynamics can be expressed in polynomial form. We also 
propose the design of decentralized control laws when such 
a certification is not possible. The approach presented here is 
parallel and scalable. Similar method can also be applicable 
to complex real world systems, such as the power system. 
While power system dynamics are non-polynomial, if those 
are transformed into polynomial forms, by introducing ad¬ 
ditional equality constraints [14], the methods developed in 
this article can be applied. Future work need to address the 
issues such as including bounds on the control effort, and 
relaxing the requirement of monotonic decrease of subsystem 
Lyapunov functions along the flow. 


Appendix 

A. Proof of Lemma [7] 

We note that since lim^ + 00 e^ = 0, V;, 


V8 G (0,minef 


3K, s.t. ef < 8 Vk> K,Vi. 


Let us assume, without any loss of generality, that 


3fo > 0 , s.t. x(t 0 ) G <x€l" 


Then, 

Vi, V,(t) = e,° + [ Vi(x)dz, Vt > t 0 
Jt 0 


(28) 


( 29 ) 


(30a) 


• 3 tf < t 0 + (e/ - e,°) /r), r) := sup V-fx) < 0 (30b) 

xe®} 


s.t. Vj(t) < e/, Vt > t} 


(30c) 


Hence we can argue that, 

Vi{t) < e,°, Vt > to,Vi 

=» 3 f 1 := maxf/, s.t ■ Vi(t) < E-,Vt > t 1 ,Vi (31) 

i 

Following similar arguments it is easy to show that, 

V{(t) < £°, Vt > t 0 , Vi 

=» Vk , 3 t k > to, s.t. Vi(t) < ef, Vt > t k , Vi (32) 
Finally combining d28l > and d32l) we observe, 

V 8 G ^0,min£° 
which concludes the proof, because of (fl2l >. 
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